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1.  Introduction 


Typically,  state  estimators  are  designed  to  blend  the  information  present  in  a  model 
of  how  the  states  progress  through  time,  called  the  process  model,  with 
measurements  of  the  system  outputs.  For  example,  consider  the  case  of  a  linear 
discrete-time  Kalman  filter  for  estimating  a  linear  time-invariant  system  with  a 
process  model  given  by: 

x*=Fxt_1+BuJt_1+wi_1.  (!) 

Where  the  noise  input  w(  is  a  white  random  sequence  and  the  deterministic  input 
vector  ut  ,  is  known.  This  model  is  used  to  form  the  state  prediction  i~+1  from  the 
last  corrected  state  estimate  xA_, .  That  is: 

xA:  =  Fi£_1+But_1.  (2) 

The  Kalman  filter  seeks  to  blend  this  prediction  with  a  measurement  of  the  current 
state,  which  is  corrupted  by  the  noise  term  vk : 

Zk  =  Hxk  +vt .  (3) 

It  accomplishes  this  task  by  computing  a  Kalman  gain  K  so  that  when  the  corrected 
state  is  calculated  via: 

xf  =  xf  +  K  (zj  -  Hxj ) ,  (4) 

the  covariance  of  the  new  estimate  error  is  minimized.1  In  some  applications 
however,  there  are  no  measurements  present,  but  there  is  more  than  one  way  to 
predict  the  state  based  on  the  previous  estimate.  Suppose  there  is  a  second  process 
model: 

X*  =Gxh  +Ent_,+  vt_,,  (5) 

with  noisy  inputs  vA  and  known  deterministic  inputs  nA  , .  The  problem  is  to  combine 
the  2  predictions  xA  /  =  FxA  ,  +  Bu  and  i~k  =  Gxt  +En/f so  that  the  combined  state 

estimate  has  the  smallest  variance.  This  situation  does  not  fit  the  format  for  which 
observability  is  defined.  At  no  point  is  the  true  state  \k  mapped  to  a  measurement. 

This  is  different  from  the  multiple  model  estimation  problem,  which  uses 
measurement  residuals  to  make  an  optimal  guess  at  which  process  model  is  correct. 
The  problem  discussed  here  is  most  similar  to  decentralized  Kalman  filtering, 
which  combines  several  independent  Kalman  filter  estimates,  but  all  of  the  Kalman 
filters  use  the  same  propagation  equations.  The  main  issue  tackled  in  decentralized 
filtering  is  dealing  with  networks  of  sensors  that  may  or  may  not  have  access  to  the 
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“central”  state  estimate  and  covariance  matrix  and  the  various  timing  and  data 
transfer  issues  associated  with  such  systems. 

This  problem  could  present  itself  in  navigation  applications.  In  the  example  used 
here,  one  process  model  comes  from  projectile  flight  dynamics,  while  another 
process  model  comes  from  inertial  measurement  unit  (IMU)  integration.  The 
primary  motivation  for  this  work  is  to  improve  on  the  methodologies  used  in  Fairfax 
and  Fresconi,2  which  addresses  this  application  in  more  detail.  Another  application 
would  be  the  blending  of  2  uncorrected  inertial  measurement  units,  which  is 
becoming  more  prevalent  as  the  size  and  cost  of  microelectromechanical  system 
(MEMS)  inertial  devices  continues  to  decrease. 

In  the  next  section,  an  estimator  is  derived  for  the  simple  linear  2-process-model 
system  described  in  the  introduction.  In  Section  3  a  simulation  example  is  given 
that  blends  a  dynamic  model  of  a  ballistic  projectile  with  an  IMU  output.  Section  4 
discusses  the  conclusions  and  future  work. 

2.  Estimator  Derivation 


The  simple  scenario  involves  process  model  1,  which  is  given  above  in  Eq.  1,  and 
has  an  initial  state  error  covariance  of: 

PA+_U  =E[et_1/e[_1>/],  (6) 

and  process  noise  that  is  normally  distributed  as: 

0,QW).  (7) 

The  error  vector  ek  is  defined  as  the  true  state  minus  the  estimated  state: 

e*=x*-V  (8) 

This  must  be  combined  with  process  model  2,  which  is  given  above  in  Eq.  5  and 
has  an  initial  state  error  covariance  of: 

=  E  [e*-l,*e*-l,g  ]  ’  (9) 

and  process  noise  that  is  normally  distributed  as: 

vh  ~  7V(0,Qv  ) .  (10) 

An  estimator  that  blends  these  2  process  models  is  hypothesized  to  have  the  fonn: 

K  =Fik,  +B«b-i  +  K(Gx;_1  +EnM  -  Fxy,  -Buw).  (11) 
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The  goal  is  then  to  find  what K  must  be  to  minimize  E[e+e*7  ]  .  The  expression  for 

the  estimator  error  can  be  detennined  by  subtracting  the  estimated  state  from  the 
true  state: 


e,  =  x,  -x. 


(12) 


=  Fe(_i  /  +w,  -K(Gx;_,  +Eni_,  - Fij_, -Buw)  ' 

The  term  inside  the  parenthesis  can  be  written  in  terms  of  the  errors  by  noting  that: 

~xk  +  xa  =  -Gx,_,  -Ent_,  -  v,_,  +Fxt_1  +But_,+w1_1  =  0  (13) 

The  new  expression  for  the  estimator  error  is: 

<  =FeU/  +wa-i  +K(Ge;_,  +vt_,  -Fey,  (14) 

Using  the  assumption  that  the  process  noise  from  the  2  systems  is  uncorrelated  leads 
to  the  simplification: 


E  wvr  =  E  vwr  =  0 . 


(15) 


A  second  assumption  is  that  noise  terms  are  uncorrelated  with  the  current  state 
errors: 


E  e*v[]  =E  e*w[]  =  0  . 


(16) 


There  are  2  situations  under  which  the  estimator  is  being  used.  The  first  situation  is 
that  the  estimator  has  been  continually  blending  the  2  process  models.  In  this  case 
the  errors  at  the  previous  step  are  identical,  and  have  the  same  covariance,  that  is: 


/rfP+  p+T  "1  _  p+ 

^  |_c*-lcjt-l  J  k—l 


(17) 


The  estimator  covariance  is  then: 


E[eyey]  =  FPyiFr+Qw 

. . .  +  (-FP;_,Fr  -  Q ,  +  FP;_,Gr ) Kr  +  K (-FP;_,Fr  - Qw  +  GP<+_,Fr ) .  (18) 

...  +  K (FP;_,Fr  +  GP;_,G T  - GP;_,Fr  - FP;_,Gr  +  Q„,  +  Qv ) kt 

The  value  of  K  is  determined  by  taking  the  derivative  of  the  trace  of  the  estimator 
covariance  and  solving  for  the  zero  slope  condition,  that  is: 


d[?r(£[eXr]) 


dN 


=  2K  (FPy,Fr  +  GP;  ,Gr  -  GPy,F7'  -  FP;_,Gr  +  Qm  +  Qv )  _  ( 1 9) 


,..-2(FP;_1Fr+Qw-FP;_1Gr) 
The  estimator  gain  is  therefore: 
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K  =  (FP;_,Fr  +Q„  -FP;_,Gr)(FP;_1Fr  +  GP;„,Gr  -  GP;_,Fr  -  FP;_,Gr  +  Q„  +  Qv )"' .  (20) 

Substituting  this  into  parts  of  Eq.  12  simplifies  the  estimator  covariance  expression 
to: 

P*+  =FPA21Fr+Q„,-K(FPJ21Fr+Q,i,-GPi21Fr)  .  (21) 

In  some  cases  the  2  process  models  are  identical,  such  as  when  blending  2  IMUs. 
This  simplification  reduces  the  gain  to  a  simple  function  of  the  process  noise 
covariance  matrices: 

K  =  QW(Q„  +  QV) (22) 

The  covariance  matrix  update  is  then: 

P*+  =  FP^F"  +  Q„  -  Qw  (Q„  +  Qv  )_1  Q„ .  (23) 


The  second  situation  is  that  the  process  models  have  been  running  independently 
from  each  other  up  until  the  current  time.  In  this  case,  their  errors  are  assumed  to 
be  uncorrelated:3 


^[ei-i.ye^]  =  =  0  . 

(24) 

The  estimator  gain  is: 

W 

n 

yd 

,/F7'  +Q11.)(FP;_1,/Fr  +GP;.ljgGr  +Q„.  +QV)  1 
\  -1 

(25) 

ii 

yd 

yc' 

+pfJ 

The  covariance  update  is: 

yd 

+ 

ii 

yd 

- 

(26) 

3.  Simulation  Example 

To  verify  the  results  in  the  above  section,  a  simulation  experiment  was  performed 
on  2  simple  process  models  with  known  statistics.  The  first  is  a  point  mass  flight 
dynamic  model  of  a  ballistic  projectile  that  is  driven  by  a  white  noise  forcing 
function.  The  next  model  is  a  simplified  inertial  measurement  model,  in  which  the 
accelerations  of  the  projectile  are  numerically  integrated  after  adding  white  noise 
to  them. 

3.1  Flight  Dynamic  Model 

The  dynamic  model  for  the  ballistic  projectile  has  the  form: 
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0 

+  1  1  +  W  . 

UJ 


(27) 


p^x2  +  z2 ACx0  Jx 


2  m 


Where  the  following  constants  are  defined  below  in  the  table: 

Table  Nomenclature  and  values  used  for  simple  dynamic  model 


Name 

Symbol 

Value  Used 

Units 

Air  Density 

P 

1.2 

kg/m3 

Cross-Sectional  Area 

A 

ti0.1552/4 

2 

m 

Axial  Force  Coefficient 

Cx  0 

0.2 

nondimensional 

Mass 

m 

40 

kg 

Gravitational 

g 

9.81 

m/s2 

Acceleration 

The  initial  conditions  for  the  projectile  state  vector  are  given  by: 


0  m 

z 

0  m 

X 

500  m/s 

z 

-500  m/s 

(28) 


The  dynamic  model  is  nonlinear.  Normally  the  equations  of  motion  from  Eq.  27 
would  be  integrated  numerically  with  a  Runge-Kutte  integration  algorithm. 
However,  to  keep  the  statistics  more  “known”  in  both  models,  a  simple  Euler 
integration  algorithm  is  used  to  fonn  the  discrete  time  model: 


4-  i 

xk 

xk- 1 

4-i 

zt 

Z  k-\ 

p  Jx2  ,  +  zt  ,  ACrn 

< 

>  =  • 

• 

r  V  K  —  1  K  —  1  XU 

- Xlr  1 

•  +  wt._1 

4 

4- 1 

K  —  1 

2m 

.4. 

A-i. 

Pa/4-1  +  4-i^  c  o  . 

zk-\  +  g 

_ 

2m 

(29) 


In  order  to  propagate  the  covariance  of  the  model,  the  state  transition  matrix  is 
formed  by  taking  the  jacobian  of  the  state-propagation  equations  with  respect  to  the 
states: 


l  o 
o  1 


F(A4-,= 


o  o 


0  0 


At 

0 

AtpX^  +  ^Cx0  Atpz^_t 

+44 

ACx0  Atpjy,  zk_x 
2w1/ii2_i  +ij_i 


0 

At 


^Cx0  /SXpxk_lzk_x 

+44 

1  ^Cx0Atp  xtl  +2ACx0Atpz/[_I 

+  4-i 


.  (30) 
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The  process  noise  covariance  matrix  was  determined  by  using  the  following 
equation,  although  the  resulting  expression  is  too  lengthy  to  display: 

Qw  =  J*  F( At  -  x)  Q  Fr  (A t-x)dx.  (31) 


The  Q,  matrix  is  given  by:4 


Qt  = 


o  o 
0  0 
0  0 
0  0 


0 

0 

erf,. 

0 


At . 


(32) 


An  example  10-second  trajectory  is  plotted  in  Fig.  1. 


Fig.  1  Example  dynamic  model  ballistic  trajectory 


An  ensemble  of  100  “true”  trajectories  was  run  with  a,,,  =  1 0  m/s2  by  integrating 
Eq.  29.  An  estimated  trajectory  was  run  by  integrating: 


X*,/  =  **-!,/  +< 


P"\/a-i  +  A-iTCjo  f 


2m 

pVa-i  +  A-iTCrl 


2  m 


A-i+g 


>A t . 


(33) 
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The  resulting  root  sum  squared  (rss)  deviation  from  the  nominal  trajectory  run  with 
zero  noise  was  compared  to  the 1(7  values  taken  from  the  diagonals  of  the  propagated 
error  covariance  matrix.  The  2  results  are  shown  in  Fig.  2.  The  figure  shows  that 
the  covariance  propagation  is  modeled  consistently  with  the  random  process  being 
simulated. 


0  5  10  0  5  10 

time  (sec)  time  (sec) 


Fig.  2  Predicted  vs.  measured  error  standard  deviation  for  the  dynamic  model 


3.2  Inertial  Measurement  Model 


The  inertial  measurement  model  consists  of  integrating  accelerometer  signals  ax  and 

a: : 


The  true  accelerometer  values  ax  and  a.  are  generated  from  the  dynamic  model 
equations: 


(35) 


Each  simulation  of  the  dynamic  model  therefore  also  produces  a  noiseless 
accelerometer  trajectory,  which  is  then  corrupted  with  noise: 
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(36) 


where  each  noise  term  vrA  and  v_Ais  normally  distributed  with  zero  mean  and 
standard  deviation  ctv  . 


Although  the  accelerometers  are  generated  using  the  dynamic  model  equations,  in 
practice  the  dynamic  model  is  usually  not  known  or  utilized  when  integrating  an 
IMU.  Using  this  modeling  approach,  the  accelerometers  act  as  deterministic  inputs: 


X 


k,g 


G(At) 


'  0 

0  " 

'  0 

0 

0 

i  M 

0 

At 

0 

KJ 

0 

0 

At 

Atg 

(37) 


G(A/)  = 


1  0 
0  1 
0  0 
0  0 


At 

0 

1 

0 


0 

At 

0 

1 


(38) 


Similar  to  the  dynamic  model  case,  the  discrete  time  process  noise  covariance  is 
calculated  from: 

Qv  =  j^G(A^-i)Q..Gr  [At-x)dx  .  (39) 

The  Q  matrix  is  given  by: 


0  0 


0  0 


o 

o 

Tv 

0 


0 

0 


(40) 


Each  estimated  trajectory  from  the  IMU  model  is  calculated  by  integrating: 


"  0 

0" 

'  0 

0 

0 

fa  1 

0 

At 

0 

IAJ 

0 

0 

At 

Atg 

(41) 


The  ensemble  standard  deviation  is  compared  to  the  standard  deviation  from  the 
propagated  covariance  when  av  =  10  m/s2.  The  results  are  shown  in  Fig.  3. 
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Fig.  3  Predicted  vs.  measured  error  standard  deviation  for  the  IMU  model 


3.3  Estimator  Performance 

The  multiple  process  model  estimator  (MPME)  from  Section  2  was  used  to 
combine  the  predicted  dynamic  model  trajectory  from  Eq.  33  with  the  predicted 
IMU  trajectory  from  Eq.  41.  Figure  4  verifies  that  the  covariance  update  given  in 
Eq.  15  is  consistent  with  the  actual  measured  error  covariance. 
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Fig.  4  Predicted  vs.  measured  error  standard  deviation  for  the  state  estimator 

Fig.  5  shows  that  for  this  example,  the  MPME  produces  a  lower  ensemble  variance 
than  either  independent  process  model.  The  simulation  example  presented  in  this 
section  was  also  run  with  different  values  of  ctw  and  a,  .  In  all  cases  the  estimator  has 

less  error  than  either  independent  process  model  on  average,  but  the  improvement 
is  more  noticeable  the  closer  the  2  noise  values  are  to  each  other.  The  MPME  was 
also  run  with  the  process  models  flipped — that  is, 

\)  =  Gif,  +  En,  +k(Fx);_1  +Bu/t_,  -Gxf,  -  En  , )  — with  identical  results.5 
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Fig.  5  Comparison  of  the  ensemble  standard  deviation  values  from  the  2  independent 
process  models  and  the  combined  state  estimate 


4.  Conclusions 


A  multiple  process  model  estimator  was  designed  to  blend  the  predictions  of  2 
separate  discrete -time  process  models  when  no  measurements  are  available  to 
obtain  an  optimal  state  estimate.  The  estimator  was  validated  with  a  simple  example 
with  known  statistics  in  which  one  process  model  based  on  projectile  flight 
dynamics  was  combined  with  a  separate  process  model  based  on  inertial 
measurement  integration.  Although  the  systems  were  greatly  simplified,  they  were 
useful  for  verifying  the  estimator,  ft  was  found  that  on  average,  the  estimator 
produces  state  estimates  with  lower  variance  than  either  independent  process 
model.  These  results  should  prove  useful  in  situations  where  the  process  is  truly 
random  and  the  statistical  distributions  of  the  forcing  functions  are  known. 

The  case  where  the  dynamic  model  is  not  actually  driven  by  white  noise  was  not 
addressed  here.  It  is  expected  that  the  marginal  improvements  gained  by  blending 
2  noisy  process  models  will  be  lost  if  either  of  the  models  includes  deterministic 
modeling  errors. 
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3.  The  -  superscript  denotes  that  the  errors  have  not  been  through  any  kind  of 
“corrector”  step. 

4.  The  A?term  comes  from  the  fact  that  Euler  integration  is  being  used  to 
propagate  the  model. 

5.  The  entire  systems  where  switched  including  process  noise  covariance 
matrices.  The  point  of  the  exercise  was  to  show  that  it  doesn’t  matter  which 
process  model  is  considered  model  1  and  which  is  model  2. 
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